• Неметрическое многомерное шкалирование
  • Как работает неметрическое многомерное шкалирование
  • Оценка качества подгонки ординации
  • Трактовка результатов ординации

Вы сможете

  • Построить диаграмму nMDS.
  • Охарактеризовать качество ординации с помощью величины стресса.
  • Сравнить результаты нескольких ординаций.

Проблема изображения многомерных данных

Облако точек в многомерном пространстве

Когда признаков много, можно представить все объекты как облако точек в многомерном пространстве.

Migration by Don McCullough on Flickr

Для изображения N объектов в идеале нужно N-1 измерений

  • 2 объекта = 1D прямая
  • 3 объекта = 2D плоскость
  • 4 объекта = 3D объем
  • N объектов = (N-1)-мерное пространство

Многомерное пространство сложно изобразить. Есть два пути:

  • Выбрать наиболее информативную проекцию (не все можно хорошо спроецировать).
  • Сохранить отношения между объектами (придется исказить расстояния).

Ежик в тумане

black shadows for a white horses / les negres ombres dels cavalls blancs by Ferran Jordà on Flickr

Неметрическое многомерное шкалирование

Nonmetric Multidimensional Scaling, nMDS (Shepard 1962, 1966, Kruskal 1964a, b)

Метод визуализации отношений между объектами в пространстве с небольшим числом измерений (обычно 2).

Исходные данные — матрица расстояний между объектами в многомерном пространстве.

nMDS ординация старается сохранить отношения между объектами.

Карта пригородов Санкт-Петербурга

Расстояния по автодорогам

Расстояния по автодорогам не совсем “евклидовы”: например, из Кронштадта до Санкт-Петербурга нельзя добраться по прямой.

Обычная проекция исходного пространства на плоскость неудачна

Карта







Проекция при помощи анализа главных координат (PCoA)

  • Карта неудачна: Стрельна, Красное Село и Гатчина сливаются вместе, хотя они не рядом.

Можно восстановить карту, сохраняя ранги расстояний между объектами

Карта





Ординация nMDS

  • Что странного в этой карте?

Важные свойства nMDS

  1. Ординация сохраняет ранг расстояний между объектами (похожие располагаются близко, непохожие — далеко; если объект А похож на B, больше чем на C, то и на ординации он, скорее всего, окажется ближе к B, чем к C).

  2. На ординации nMDS имеет смысл только взаиморасположение объектов. Облако точек в осях MDS можно вращать, перемещать, зеркально отражать. Суть ординации от этого не изменится.

  3. Численные значения координат точек в ординации лишены смысла (их вообще можно не приводить на итоговой ординации).

Как работает nMDS

Наблюдаемые различия и ожидаемые расстояния

Взаиморасположение точек на плоскости подобно взаиморасположению точек в многомерном пространстве признаков, но значения не полностью совпадают.

Диаграмма Шепарда

Диаграмма Шепарда — график зависимости расстояний между точками на ординации и соответствующих им значений коэффициентов различия в исходном пространстве.


  • По оси X — значения коэффициентов различий \(D_{h,i}\), ранжированные в порядке возрастания
  • По оси Y — значения расстояний на плоскости ординации — \(d_{h,i}\)

В идеале, ранги расстояний между точками \(d_{h,i}\) на хорошей ординации nMDS должны совпадать с рангами коэффициентов различия \(D_{h,i}\) в исходном многомерном пространстве.

Интерпретация диаграммы Шепарда

Форма и расположение облака точек на диаграмме Шепарда позволяет судить о качестве ординации.


  • A — ординация объясняет большую часть исходной изменчивости
  • В — ординация объясняет небольшую часть изменчивости, но относительное расположение точек сохранено
  • С — относительное расположение точек плохо соответствует исходному

Диаграмма Шепарда

В идеале, большим расстояниям в исходном пространстве должны соответствовать большие расстояния на ординации, а меньшим — меньшие (т.е. должны сохраняться ранги расстояний).


Красная линия — предсказанные монотонной регрессией значения расстояний на ординации — монотонно возрастает.

Монотонная регрессия

Монотонная регрессия минимизирует сумму квадратов отклонений значений координат на ординации от значений в исходном многомерном пространстве, сохраняя монотонные отношения между ними.

Алгоритм подбора монотонной регрессии называется Pool-Adjacent-Violator-Algorithm (PAVA)

Стресс — мера качества ординации

Стресс (Stress) показывает, насколько соответствуют друг другу взаиморасположение точек на плоскости ординации и в исходном многомерном пространстве признаков.

Стресс 1 (Stress(1), Kruskal, 1964a) — это корень из суммы квадратов остатков монотонной регрессии, отнесенной к сумме квадратов всех коэффициентов различия в исходной матрице.

\[ Stress(1) = \sqrt{\frac{\sum_{h < i}{(d_{h,i} - \hat d_{h,i})^2}}{\sum_{h < i}{d_{{h,i}}^2}}}\]

  • \(d\) — наблюдаемое значение коэффициента различия между двумя точками
  • \(\hat d\) — значение, предсказанное монотонной регрессией

Хорошо, когда \(Stress(1) < 0.2\).

Алгоритм MDS (для двумерного случая)

  1. Выбираем a priori число измерений k для финального решения (число осей пространства ординации).
  2. Вычисляем матрицу коэффициентов различия между объектами.
  3. Выбираем начальное расмещение объектов в пространстве k осей (случайно или при помощи PCoA).
  4. Вычисляем матрицу расстояний между точками в пространстве ординации.
  5. Вычисляем стресс.
  6. Сдвигаем точки, чтобы минимизировать стресс (метод steepest descent).
  7. Повторяем шаги 3-6 несколько раз, чтобы избежать локальных минимумов стресса.
  8. Обычно финальную ординацию поворачивают при помощи PCA, чтобы вдоль оси X было максимальное варьирование.

Ежик в тумане

Ограничения nMDS

  • nMDS — метод визуализации данных.
  • Нужно осознанно выбирать способ трансформации данных и коэффициент различия исходя из того, какие именно свойства отношений между объектами должны быть визуализированы.
  • Число измерений, выбранное для решения, может влиять на результат.

Пример

Симбионты мидий

Открываем данные

dat <- read.delim("data/Krapivin_2017_Medulis-symb.tab", skip = 36, 
                  stringsAsFactors = TRUE)
str(dat)
## 'data.frame':    548 obs. of  23 variables:
##  $ Event                 : Factor w/ 3 levels "Gorelyi_Isl",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ Date.Time             : Factor w/ 1 level "2013-07": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Ord.No                : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ Area                  : Factor w/ 3 levels "B Solovkecki Isl",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Site                  : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ Zone                  : Factor w/ 3 levels "down","middle",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Comment               : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ M..edulis.L.shell..mm.: num  31.2 41.1 46.4 40 34.5 47.6 34.3 49.6 49.8 30 ...
##  $ M..edulis.age..a.     : num  1 2 2 3 4 3 2 4 3 1 ...
##  $ Urastoma.spp.....     : int  26 16 19 11 14 16 10 4 11 0 ...
##  $ Renicola.spp.....     : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Himasthla.spp.....    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Metacercaria....      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Gymnophallus.spp..... : int  0 0 0 0 0 0 0 0 4 0 ...
##  $ Bucephalidae....      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Alg....               : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Nematoda....          : int  1 1 0 0 0 0 0 0 0 0 ...
##  $ Microsetella....      : int  0 0 0 0 0 2 0 0 0 0 ...
##  $ Copepoda....          : int  0 0 0 0 0 1 0 0 0 0 ...
##  $ Chironomidae....      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Halacaridae....       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Jaera.spp.....        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Ostrac....            : int  0 0 0 0 0 0 0 0 0 0 ...

Приводим в порядок названия переменных

colnames(dat) # Было
##  [1] "Event"                  "Date.Time"              "Ord.No"                
##  [4] "Area"                   "Site"                   "Zone"                  
##  [7] "Comment"                "M..edulis.L.shell..mm." "M..edulis.age..a."     
## [10] "Urastoma.spp....."      "Renicola.spp....."      "Himasthla.spp....."    
## [13] "Metacercaria...."       "Gymnophallus.spp....."  "Bucephalidae...."      
## [16] "Alg...."                "Nematoda...."           "Microsetella...."      
## [19] "Copepoda...."           "Chironomidae...."       "Halacaridae...."       
## [22] "Jaera.spp....."         "Ostrac...."
colnames(dat) <- gsub("[.]", replacement = "", colnames(dat))
dat <- rename(dat, L = "MedulisLshellmm", Age = "Medulisagea")

colnames(dat) # Стало
##  [1] "Event"           "DateTime"        "OrdNo"           "Area"           
##  [5] "Site"            "Zone"            "Comment"         "L"              
##  [9] "Age"             "Urastomaspp"     "Renicolaspp"     "Himasthlaspp"   
## [13] "Metacercaria"    "Gymnophallusspp" "Bucephalidae"    "Alg"            
## [17] "Nematoda"        "Microsetella"    "Copepoda"        "Chironomidae"   
## [21] "Halacaridae"     "Jaeraspp"        "Ostrac"

Приводим в порядок данные

# Делаем сайт фактором
dat$Site <- factor(dat$Site, levels = c(2, 3, 1), labels = c("G","L","S"))

# Сливаем редкие виды в одну категорию
f_remove <- c("Nematoda", "Microsetella", "Copepoda", "Chironomidae", 
              "Halacaridae", "Jaeraspp", "Ostrac")
dat$Other <- rowSums(dat[, f_remove])

# Суммарная численность симбионтов
f_sp <- c("Urastomaspp", "Renicolaspp", "Himasthlaspp", "Metacercaria", 
          "Gymnophallusspp", "Alg", "Other")
dat$Total  <- rowSums(dat[, f_sp])

# Данные для анализа
# Только мидии с симбионтами и возрастом от 3 до 8 лет 
dat <- dat[dat$Total != 0 & dat$Age %in% 3:8, ]
spec <- dat[, f_sp]                         # виды-симбионты
env <- dat[, c("Zone", "Site", "L", "Age")] # свойства мидий-хозяев

Ординация сообществ симбионтов в мидиях

Задание 1

Постройте ординацию мидий по обилию разных видов-симбионтов с использованием коэффициента различия Брея-Куртиса. Следите, чтобы алгоритму удалось найти финальное решение. Если необходимо, настройте вызов metaMDS()

Вычислите стресс для получившейся ординации.

Нарисуйте простейший график.

ord_mus <- 

Решение

ord_mus <- metaMDS(spec, dist = "bray")

По-умолчанию metaMDS() пытается подобрать оптимальную трансформацию, подходящую для анализа сообществ. В данном случае — извлекает корень и делает двойную висконсинскую стандартизацию.

Square root transformation
Wisconsin double standardization
Run 0 stress 0.1319177 

Такая модель не сходится.

*** No convergence -- monoMDS stopping criteria:
    12: stress ratio > sratmax
     8: scale factor of the gradient < sfgrmin

Решение

После отключения автоматического подбора трансформации данных удается найти решение.

ord_mus <- metaMDS(spec, dist = "bray", autotransform = FALSE)
## Run 0 stress 0.141 
## Run 1 stress 0.145 
## Run 2 stress 0.144 
## Run 3 stress 0.147 
## Run 4 stress 0.142 
## Run 5 stress 0.141 
## ... New best solution
## ... Procrustes: rmse 0.00132  max resid 0.0123 
## Run 6 stress 0.163 
## Run 7 stress 0.153 
## Run 8 stress 0.141 
## ... Procrustes: rmse 0.00181  max resid 0.0326 
## Run 9 stress 0.141 
## ... Procrustes: rmse 0.000133  max resid 0.002 
## ... Similar to previous best
## Run 10 stress 0.141 
## ... Procrustes: rmse 0.00339  max resid 0.0334 
## Run 11 stress 0.155 
## Run 12 stress 0.141 
## ... Procrustes: rmse 0.00121  max resid 0.0119 
## Run 13 stress 0.141 
## ... Procrustes: rmse 0.0014  max resid 0.0119 
## Run 14 stress 0.156 
## Run 15 stress 0.156 
## Run 16 stress 0.141 
## ... New best solution
## ... Procrustes: rmse 0.00062  max resid 0.0093 
## ... Similar to previous best
## Run 17 stress 0.141 
## ... Procrustes: rmse 0.00365  max resid 0.0337 
## Run 18 stress 0.146 
## Run 19 stress 0.141 
## ... Procrustes: rmse 0.00235  max resid 0.0335 
## Run 20 stress 0.146 
## *** Solution reached

Решение

Хорошее решение, но график нужно украсить.

ord_mus$stress
## [1] 0.141
ordiplot(ord_mus)

Украшенный график

# Палитры
pal_col <- c("red", "green", "steelblue")
pal_sh <- c(1, 2, 0)

# Украшенный график
ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])

Украшенный график с центроидами видов

ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
text(ord_mus, display = "spec", cex = 0.9, col = "grey20")

Визуализация ординации в ggplot2

Задание 2

Используя данные, приведенные ниже, постройте график ординации при помощи ggplot2. Покажите точками — мидий, текстом — центроиды для видов-симбионтов, цветом — зону литорали, формой — сайт, размером маркеров — размер мидий.

# Координаты точек (мидий)
ord_mus_pt <- data.frame(env, scores(ord_mus, display = "sites"))
head(ord_mus_pt, 2)
# Координаты центроидов переменных (видов-симбионтов)
ord_mus_sp <- data.frame(scores(ord_mus, display = "species"))
ord_mus_sp$Species <- rownames(ord_mus_sp)
head(ord_mus_sp, 2)

Решение: График с точками

gg_ord_mus <- ggplot() +
  geom_point(data = ord_mus_pt, 
             aes(x = NMDS1, y = NMDS2, colour = Zone, 
                 shape = Site, size = L), alpha = 0.5)
gg_ord_mus

Решение: График с центроидами видов

gg_ord_mus_sp <- gg_ord_mus +
  geom_text(data = ord_mus_sp, 
            aes(x = NMDS1, y = NMDS2, 
                label = Species))
gg_ord_mus_sp

Интерпретация ординации: envfit

Анализ связи с внешними переменными c помощью функции envfit()

Функция envfit() строит регрессионную модель зависимости переменной среды от координат точек на плоскости ординации.

Модель вот такого вида:

\(y_i = \beta_1 \text{NMDS1}_i + \beta_2 \text{NMDS2}_i + \varepsilon_i\)

Значимость зависимости оценивается при помощи пермутационного теста (это позволяет работать в ситуации, когда не полностью выполнены условия применимости).

Такой вызов функции позволит подобрать сразу несколько моделей и изобразить их на графике ординации:

ef <- envfit(ord_mus, env[, c("Zone", "Site", "L", "Age")])

Интерпретация результатов envfit() для непрерывных переменных

ef$vectors
##      NMDS1  NMDS2   r2 Pr(>r)    
## L   -1.000 -0.022 0.19  0.001 ***
## Age  0.777 -0.629 0.04  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


  • Колонки MDS1 и MDS2 содержат информацию о направлении вектора непрерывной переменной. Длина вектора пропорциональна корню из коэффициента детерминации для модели.
  • r2 — \(R^2\) из регресcионного анализа.
  • Pr(>r) — пермутационная оценка статистической значимости.

Интерпретация результатов envfit() для дискретных переменных

ef$factors
## Centroids:
##            NMDS1 NMDS2
## Zonedown   -0.65 -0.09
## Zonemiddle -0.14 -0.05
## Zoneup      0.83  0.14
## SiteG       0.19  0.17
## SiteL       0.23 -0.41
## SiteS      -0.42  0.23
## 
## Goodness of fit:
##        r2 Pr(>r)    
## Zone 0.37  0.001 ***
## Site 0.16  0.001 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999


  • Колонки MDS1 и MDS2 содержат координаты центроидов для дискретных факторов (т.е. положение центра облака точек, относящихся к категории фактора).
  • r2 — \(R^2\) из регресcионного анализа.
  • Pr(>r) — пермутационная оценка статистической значимости.






График с векторами и центроидами, найденными envfit()

ordiplot(ord_mus, type = "n")
points(ord_mus, col = pal_col[env$Zone], pch = pal_sh[env$Site])
plot(ef)

Для ggplot-графика удобно использовать вспомогательный пакет

# install.packages("devtools")
# devtools::install_github("gavinsimpson/ggvegan")
library(ggvegan)
ord_mus_ef <- fortify(ef)
ord_mus_ef

ggplot2 версия графика envfit()

gg_ord_mus +
  geom_segment(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
               aes(x = 0, xend = NMDS1, y = 0, yend = NMDS2),
               arrow = arrow(length = unit(0.25, "cm"))) +
  geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Vector", ],
            aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1)) +
  geom_text(data = ord_mus_ef[ord_mus_ef$Type == "Centroid", ],
            aes(x = NMDS1, y = NMDS2, label = Label, hjust = 1.1, vjust = 1))

Ограничения envfit()

  • Осторожно! envfit() основан на предположении о линейном характере связи между внешней переменной и координатами на ординации. Это не всегда так.

Интерпретация ординации: ordisurf

Анализ связи с внешними переменными c помощью ordisurf()

ordisurf() применяется, когда зависимость нелинейна.

В основе работы функции — общая аддитивная модель (General Additive Model, GAM).

GAM позволяют подобрать нелинейные зависимости, даже если заранее неизвестна их форма.

ordisurf() использует пакет mgcv для подбора GAM на основе метода тонких пластин (thin plate spline).

График с поверхностью, найденной ordisurf()

Изолинии на графике ординации показывают предсказанное значение зависимой переменной для разных участков ординации.

par(mfrow = c(1, 2))
os_L <- ordisurf(ord_mus, env$L, method = "REML")
os_Age <- ordisurf(ord_mus, env$Age, method = "REML")
par(mfrow = c(1, 1))

Интерпретация результатов ordisurf()

Длина мидий линейно меняется на ординации (эффективное число степеней свободы для сплайна \(edf=1.96\)). Зависимость довольно хорошо выражена (объясненная девианса 18.7%).

summary(os_L)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   43.410      0.484    89.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df    F p-value    
## s(x1,x2) 1.96      9 9.83  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.183   Deviance explained = 18.7%
## -REML = 1457.6  Scale est. = 92.442    n = 395

Интерпретация результатов ordisurf()

Возраст мидий меняется нелинейно вдоль ординационных осей (эффективное число степеней свободы для сплайна \(edf=4.99\)). Эта зависимость довольно слабая (объясненная девианса 8.97%).

summary(os_Age)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## y ~ s(x1, x2, k = 10, bs = "tp", fx = FALSE)
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   4.8861     0.0714    68.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##           edf Ref.df    F    p-value    
## s(x1,x2) 5.01      9 3.73 0.00000097 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.0785   Deviance explained = 9.02%
## -REML =  704.9  Scale est. = 2.0111    n = 395

Для ggplot-графика понадобится добыть данные о контурах

fortify_ordisurf <- function(model) {
  # Fortifies an object of class ordisurf to produce 
  # a dataframe of contour coordinates 
  # for subsequent plotting in ggplot.
  xy <- expand.grid(x = model$grid$x, y = model$grid$y)
  xyz <- cbind(xy, c(model$grid$z))
  names(xyz) <- c("x", "y", "z")
  return(na.omit(xyz))
}

ord_mus_os <- fortify_ordisurf(os_Age)
head(ord_mus_os, 4)

ggplot2 версия графика с поверхностью, найденной ordisurf()

ggplot(data = ord_mus_os, aes(x = x, y = y, z = z)) +
  stat_contour(aes(colour = ..level..)) +
  labs(x = "NMDS1", y = "NMDS2", colour = "Age")

Ограничения ordisurf()

  • Осторожно! Для использования ordisurf() важно, чтобы в ординации участвовало много объектов и они располагались на плоскости ординации более-менее равномерно. Пустоты и отскоки снижают качество предсказаний.

Финальный график

Take-home messages

  • nMDS — один из методов ординации в пространстве сниженной размерности. При ординации этим методом сохраняются ранги расстояний между объектами.
  • Мера качества ординации называется “стресс”.
  • Значения координат не имеют особенного смысла. Имеет значение только взаиморасположение точек.
  • Результат nMDS зависит от выбора меры различия.
  • Существует несколько способов визуализировать на ординации значения внешних переменных. При выборе между envfit() и ordisurf() обращайте внимание на форму связи.

Литература

  • Borcard, D., Gillet, F., Legendre, P., 2011. Numerical ecology with R. Springer.
  • Legendre, P., Legendre, L., 2012. Numerical ecology. Elsevier.
  • Oksanen, J., 2015. Multivariate analysis of ecological communities in R: vegan tutorial. http://www.cc.oulu.fi/~jarioksa/opetus/metodi/vegantutor.pdf
  • Zuur, A. F., Ieno, E. N., Smith, G. M. Analysing Ecological Data. Springer 2007

Самостоятельная работа

Задание 3

Во всех примерах необходимо:

  1. Разобраться с данными.
  2. Построить ординацию объектов (описаний, проб и т.п.).
  3. Визуализировать связь между полученной ординацией и параметрами среды.
  4. Сделать выводы о наиболее важных факторах.

Источники данных

  1. Растительные сообщества во франции La Mafragh (данные из работы de Belair et al. 1987, данные mafragh, пакет ade4).
  2. Деревья на острове Barro Colorado (данные из работы Condit et al. (2002), данные BCI, пакет vegan).